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Abstract 

The method to simulate the rescattering and energy loss of hard partons in ultrarel- 
, ativistic heavy ion collisions has been developed. The model is a fast Monte-Carlo tool 

I introduced to modify a standard PYTHIA jet event. The full heavy ion event is obtained 

■ as a superposition of a soft hydro-type state and hard multi-jets. The model is applied to 

analysis of the jet quenching pattern at RHIC. 
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1 Introduction 

One of the important tools for studying the properties of quark-gluon plasma (QGP) in ultra- 
relativistic heavy ion collisions is the analysis of a QCD jet production. The medium-induced 
energy loss of energetic partons, "jet quenching", should be very different in the cold nu- 
clear matter and QGP, resulting in many observable phenomena [1]. Recent RHIC data on 
high-pT particle production [2] (the suppression of hadron spectra and azimuthal back-to-back 
two-particle correlations, strong elliptic flow) are in agreement with the jet quenching hypothe- 
sis [3]. At LHC, a new regime of heavy ion physics will be reached at ^snn = 5.5A TeV where 
hard and semi-hard particle production can stand out against the underlying soft events. The 
initial gluon densities in Pb+Pb reactions at LHC are expected to be significantly higher than 
those at RHIC, implying a stronger partonic energy loss, observable in various new channels [4]. 

In the most of available Monte- Carlo heavy ion event generators the medium-induced par- 
tonic rescattering and energy loss are either ignored or implemented insufficiently [5]. Thus, 
in order to analyze RHIC data on high-p^- hadron production, test the sensitivity of LHC 
observables to the QGP formation, and study the corresponding experimental capabilities of 
detectors, the development of adequate and fast Monte-Carlo models for jet quenching simula- 
tion is necessary. 

In this paper we present the model of jet quenching and discuss its validation basing on 
RHIC data for high-p^ hadron spectra. In Sect. 2 we give the physics frameworks of the model. 
Section 3 describes the event-by-event simulation procedure. In Sect. 4 the generalization of 
the model to the case of "full" heavy ion event (superposition of soft hydro-type state and hard 
multi-jets) is fulfilled. In Sect. 5 the efficiency of the model is demonstrated by means of the 
numerical analysis of hadron spectra at RHIC. 



2 Physics frameworks of the model 

The detailed description of physics frameworks of the developed model can be found in a number 
of our previous papers [6, 7, 8, 9, 10]. Our approach bases on an accumulating energy loss, 
the gluon radiation being associated with each parton scattering in the expanding medium and 
includes the interference effect using the modified radiation spectrum dE/dl as a function of 
decreasing temperature T. The basic kinetic integral equation for the energy loss as a 
function of initial energy E and path length L has the form 

where I is the current transverse coordinate of a parton, dP/dl is the scattering probability 
density, dE/dl is the energy loss per unit length, A = l/{ap) is in-medium mean free path. 
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p oc is the medium density at the temperature T, a is the integral cross section for parton 
interaction in the medium. 

The coUisional energy loss due to elastic scattering with high-momentum transfer has been 
originally estimated by Bjorken in [11], and recalculated later in [12] taking also into account 
the low-momentum transfer loss resulting mainly from the interactions with plasma collective 
modes. Since the latter process does not contribute much to the total coUisional loss in com- 
parison with high- momentum scattering (due to absence of large factor ~ In (£'///£>) where 
is the Debye screening mass) and, in numerical estimates it can be effectively "absorbed" by 
means of redefinition of minimum momentum transfer t^am ~ A*d , we used the coUisional part 
associated with high- momentum transfer only [7], 

di?^"' 1 , da 



dl ATXa 

and the dominant contribution to the differential cross section 

da ^ 2T,al{t) Utt 

dt E^-my ' (33 - 2Nf) In (i/A^ci.) ^ ^ 

for scattering of a hard parton with energy E and mass rrip off the "thermal" parton with 
energy (or effective mass) mo ~ 37 <S E. Here C = 9/4, 1,4/9 for gg, gq and qq scatterings 
respectively, as is the QCD running coupling constant for Nf active quark flavors, and Aqcd 
is the QCD scale parameter which is of the order of the critical temperature, Aqcd ^ Tc ^ 
200 MeV. The integrated cross section a is regularized by the Debye screening mass squared 
n\){T) ~ 47rQ;sT^(l + Nf/6). The maximum momentum transfer t^ax — [s — {rrip + mo)^][s — 
{rrip — moY]/s where s — 2moE + ml + rrip. 

There are several calculations of the inclusive energy distribution of medium-induced gluon 
radiation using Feyman multiple scattering diagrams. The relation between these approaches 
and their basic parameters has been discussed in detail in the recent writeup of the working 
group "Jet Physics" for the CERN Yellow Report [4]. We restrict ourselves to using EDMS 
formalism [13]. In the EDMS frameworks, the strength of multiple scattering is characterized 
by the transport coefficient q — /i"^ / Xg {Xg is the gluon mean free path), which is related to the 
elastic scattering cross section a (3). In our simulations this strength is in fact regulated mainly 
by the initial QGP temperature Tq. Then the energy spectrum of coherent medium- induced 
gluon radiation and the corresponding dominant part of radiative energy loss of massless parton 
have the form [13]: 



E 



dl ttL 



duj 



y + 



In |cos(a;iri)| , (4) 



u;,^^lz{l-y+^^y^]Rln^ with R ^ J^^^ , (5) 




2 



where Ti = L/{2Xg), y — uj/E is the fraction of the hard parton energy carried away by the 
radiated gluon, and Cr = 4/3 is the quark color factor. A similar expression for the gluon jet 
can be obtained by setting Cr — 3 and proper by changing the factor in the square brackets in 
(4), see ref. [13]. The integration (4) is carried out over all energies from (^min = -^lpm = l^'o^g, 
the minimum radiated gluon energy in the coherent LPM regime, up to initial parton energy 
E. Note that we do not consider here possible effects of double parton scattering [14, 15] and 
thermal gluon absorption [16], which can be included in the model in the future. 

The simplest generalization of the formula for a heavy quark of mass ruq can be done by 
using the "dead-cone" approximation [17]: 

dE 1 dE r-(AX'^(^V'^ 

dldw^""'""' ~ (1 + (/3o;)3/2)2 dldw^"''='' ^ " Ud/ ■ 

One should mention the more recent developments on heavy quark energy loss calculations 
available in the literature [18, 19, 20], which can be also considered as further model improve- 
ments. 

The medium is treated as a boost-invariant longitudinally expanding quark-gluon fluid, and 
partons as being produced on a hyper-surface of equal proper times r [21]. In order to simplify 
numerical calculations we omit here the transverse expansion and viscosity of the fluid using 
the well-known scaling solution obtained by Bjorken [21] for a temperature and density of QGP 
at T > Tc ~ 200 MeV: 

e{ry/' ^ eor^/\ T{Tyl^ ^T,tI'\ p{t)t ^ p^r,. (7) 

The internal model parameters are the initial conditions for the QGP formation expected for 
central Au-|-Au (Pb-|-Pb) collisions at RHIC (LHC): tq, Tq and Nj. For non-central coUisions 
and for other beam atomic numbers we suggest the proportionality of the initial energy density 
£o to the ratio of nuclear overlap function and effective transverse area of nuclear overlapping [7] . 

Note that using other scenarios of QGP space-time evolution for the Monte-Carlo implemen- 
tation of the model is also envisaged. In fact, the influence of the transverse flow, as well as that 
of the mixed phase at T = Tc, on the intensity of jet rescattering (which is a strongly increasing 
function of T) has been found to be inessential for high initial temperatures Tc- On the 

contrary, the presence of QGP viscosity slows down the cooling rate, that implies a jet parton 
spending more time in the hottest regions of the medium. As a result the rescattering intensity 
increases, i.e., in fact an effective temperature of the medium gets higher as compared with the 
perfect QGP case. We also do not take into account here the probability of jet rescattering 
in nuclear matter, because the intensity of this process and corresponding contribution to the 
total energy loss are not signiflcant due to much smaller energy density in a "cold" nucleus. 

Another important element of the model is the angular spectrum of in-medium gluon ra- 
diation. Since the detailed calculation of the angular spectrum of emitted gluons is rather 
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sophisticated and model-dependent [6, 13, 15, 22, 23, 24], the simple parameterization of gluon 
angular distribution over the emission angle 6 was used: 

^ocsm^expj^ ^j, (8) 

where 9q ~ 5° is the typical angle of the coherent gluon radiation as estimated in [6]. Other 
parameterizations are also envisaged. 

3 Simulation procedure 

The model has been constructed as the fast Montc-Carlo event generator PYQUEN (PYthia 
QUENched), and the corresponding Fortran routine PYQUEN is available via Internet [25]. 
The routine is implemented as a modification of the standard PYTH1A_6.2.* jet event [26]. 
The following cvent-by-event Monte-Carlo simulation procedure is applied. 

• Generation of the initial parton spectra with PYTHIA (fragmentation off). 

• Generation of the jet production vertex at the impact parameter h according to the distribu- 
tion 

where ri^2{b,r,ilj) are the distances between the nucleus centers and the jet production vertex 
V{r cos ipjT sin ip); rmax{b,ip) < Ra is the maximum possible transverse distance r from the 
nuclear collision axis to V; Ra is the radius of the nucleus A; T^ir) = AJ pA{r,z)dz is the 
nuclear thickness function with nucleon density distribution pA{r,z); T^Aib) is the nuclear 
overlap function (see ref. [7] for detailed nuclear geometry explanations). 

• Calculation of scattering cross section a = J dt da/dt (3). 

• Generation of the displacement between i-th and {i + l)-th scatterings, k = (r^+i — Ti): 



dP 
dli 



H 

A-i(Ti+i)exp(- J X-\n + s)ds) , X-\t) ^ a{T)p{T) , (10) 



and calculation of the corresponding transverse distance, liPx/E. 

• Reducing the parton energy by coUisional and radiative loss per each i-th scattering: 

A-Etot,i = A£'col,i + A£^rad,i , (H) 

where the coUisional part is calculated in the high-momentum transfer approximation (3), 

AE^oU^-^, (12) 
'° 2mo ^ ^ 
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and the energy of a radiated gluon, uji — AE'j-ad.i, is generated according to (4) and (6) 



dl . 2as{nl)XC, 



R 



J mg=0 r 

duj ttLuj 



1 ^ y" 



d/ 1 dl 

InlcosKrOI , ^U,^o = + (^^)3/.)2 ^k-o -(13) 



• Calculation of the parton transverse momentum kick due to elastic scattering i: 

Ml, = {E- ^ f -{p--J^-^f- ml. (14) 
^ 2moi p2moi 2p' ^ ^ ' 

• Formation of the additional (in-medium emitted) gluon with the energy cui and the emission 
angle 9, relative to the parent parton determined according to the parameterization (8). 

• Halting the rescattering ii (1) the parton escapes the dense zone, oi (2) QGP cools down to 
Tc — 200 MeV, or (3) the parton loses so much energy that its Pt{t) drops below 2T(r). 

• At the end of each event, adding new (in-medium emitted) gluons to the PYTHIA parton 
list and rearrangement of partons to update string formation. 

• Formation of the final state particles by PYTHIA (fragmentation ori) . 

4 Extension of the model to simulate full heavy ion event 

The full heavy ion event is simulated as a superposition of soft hydro-type state and hard 
multi-jets. The simple approximation [2 7, 28] of hadronic hquid at "freeze-out" stage has been 
used to treat soft part of the event giving the final hadron spectrum as a superposition of a 
thermal distribution and a collective flow [29, 30, 31]. 

1. The 4-momentum p* of a hadron of mass m was generated at random in the rest frame of a 
liquid element in accordance with the isotropic Boltzmann distribution 



f{E*) oc £;V^*2 _ ^2exp {-E*/Tf), - 1< cosr < 1, < 0* < 27r , (15) 



where E* = \/p*'^ + m? is the energy of the hadron, and the polar angle 6* and the azimuthal 
angle 0* specify the direction of its motion in the rest frame of the liquid element. 
2. The spatial position of a liquid element and its local 4- velocity were generated at random 
in accordance with phase space and the character of motion of the fluid: 

/(r) = 2r/R) (0 < r < it!/), /(r^) oc exp [-(77 - y-ax)2/2(yma.)2j ^ q < $ < 27r, 

Ur — sinhKji^'^"" , ^ =, Ut = Jl + ui coshr?, Uz = Jl + ui sinh ri , (16) 

^RMRfib = 0) ^ ' ^ ' 

where -R/ is the final transverse radius of the system in a given direction. Freeze-out parameters 
of the model are kinetic freeze-out temperature Tf and maximum longitudinal, Yf^'^^, and 
transverse, Y^^^, collective flow rapidities. 
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3. Further, boost of the hadron 4-momentum in the cm. frame of the event was calculated: 



Px — P* sin 9* cos 4>* + Ur cos <l> 
Py = p* sin 9* sin 0* + Ur sin $ 



E* + 
E* + 



Ut + 1 
Ut+1 



Pz = p* COS 9* + Uz E* + 
E = E*ut + {uY'): 



{u^p*"-) 

Ut + 1 



where 



(uY^) = UrP* sin 9* cos ($ - (p*) + UzP* cos 9* 



(17) 



(18) 



Anisotropic flow is introduced here under simple assumption that the spatial ellipticity of 
"freeze-out" region, e= [y^ — x^) / iy^ + x^), is directly related to the ellipticity of the system 
formed in the region of the initial overlap of nuclei, eg = 6/2i?^. This "scaling" enables one to 
avoid introducing additional parameters and, at the same time, leads to an azimuthal anisotropy 
of generated particles due to dependence of transverse radius Rf{h) on the angle $ [28]: 



i?/(6) = i?/(6 = 0) min{y 1 - eo sin"^ $ + eo cos$, \Jl-el sin^$-eo cos$} 



(19) 



Obtained in such a way azimuthal distribution of particles is described well by the elliptic form 
for the domain of reasonable impact parameter values. 

The mean total particle multiphcity in central Au+Au (Pb+Pb) coUisions at RHIC (LHC) 
is the input parameter of the model (instead of it!/(6 = 0) we put Ra here for simplicity), the 
total multiplicity for other centralities and atomic numbers being assumed to be proportional 
to the number of nucleons-participants. We also set the Poisson multiplicity distribution and 
the following particle ratios: 

TT^ : : = 24 : 6 : 1, tt^ : 7r° = 2 : 1, : = I : I, p : n = 1 : 1 . 

The hard part of the event includes PYTHIA/PYQUEN hadronic jets generated according 
to the binomial distribution. The mean number of jets produced in AA events at a given h is 
proportional to the number of binary nucleon-nucleon sub-collisions and determined as 



NfAib, V~s) = TAA{h) j dpi j dy- 



da^?,'^ipT, yi) 



dp^dy 



(20) 



where da^^'^ipx, y/s)/dp'ldy is the cross section of corresponding hard process in pp collisions 
(at the same c.m.s. energy, ^/s, of colliding beams) with the minimum transverse momentum 
transfer p™™. The latter is another input parameter of the model. In the frameworks of our 
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approximation, partons produced in (semi)hard processes with the momentum transfer less 
than p™™ are considered as being "thermahzed" , so their hadronization products are included 
in a soft part of the event "automatically" . 

Note that we can expect some adequate results only for central and semi-central collisions, 
but not for very peripheral collisions {b ~ 2Ra) where the hydro-type description is not applica- 
ble. Besides, the very forward rapidity region (where other dynamical effects can be important) 
is beyond our treatment here. 

The extended in such a way jet quenching model has been constructed as the fast Monte- 
Carlo event generator, and the corresponding Fortran code is also available via Internet [32]. 

Let us remind in the end of this section, that ideologically our approximation is similar to 
the model of Hirano and Nara [33] . The difference is that we concentrate here on the detailed 
simulation of the parton multiple scattering in a QCD-medium (the scattering-by-scattering 
generation of parton path length and energy loss in an expanding QGP, taking into account 
the coUisional loss, Lund string fragmentation model both for hard partons and in-medium 
emitted gluons, etc.), while the treatment of the hydrodynamic part in [33] is much more 
detailed than in our simple (and therefore fast) simulation procedure. 

5 Validation of the model at RHIC, = 200 A GeV 

In order to demonstrate the efficiency of the model, the jet quenching pattern in Au-|-Au col- 
lisions at RHIC was considered. The comparison of calculated and experimentally measured 
pseudorapidity r] and transverse momentum px spectra of hadrons together with their depen- 
dence of event centrality allows the optimization of the model and specification of main model 
parameters. 

The PHOBOS data on ?7-spectra of charged hadrons [34] have been analyzed to fix the 
particle density in the mid-rapidity region and the maximum longitudinal fiow rapidity, Yf^^^ — 
3.5. For the calculation of (multi)jet production cross section, we used the factor K = 2 taking 
into account higher order corrections of perturbative QCD. The rest of the model parameters 
have been obtained by fitting PHENIX data on p^-spectra of neutral pions [35] : the kinetic 
freeze-out temperature Tf = 100 MeV, maximum transverse fiow rapidity K^"^^ = 1.25 and 
minimum transverse momentum transfer of "non-thermalized" hard process p™"' = 2.8 GeV/c. 
It was found that the nuclear modification of the hardest domain of pr-spectrum (> 5 GeV/c) 
is determined in our case only by the intensity of the medium-induced parton rescattering. This 
fact allows us to extract from the data initial conditions of the QGP formation independently 
on other input parameters: the initial temperature Tq = 500 MeV, the formation time Tq = 0.4 
fm/c and the number of active quark fiavours Nj = 2. We will see below, that setting model 
parameters as it was described above, makes it possible to reproduce the main features of jet 
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quenching pattern at RHIC: the p^-dependence of the nuclear modification factor Raa and 
two-particle azimuthal correlation function C(A</?). 

Figure 1 shows jy-distribution of charged hadrons in Au+Au collisions for different centrality 
sets. The good fit of PHOBOS data [34] is achieved excepting very forward rapidities. The px- 
distributions of 7r°-mesons obtained at PHENIX [35] is also well reproduced by our calculations, 
even for relatively peripheral collisions (Figure 2). 

Figure 3 shows the nuclear modification factor Raa for neutral pions, which is defined as 
a ratio of particle yields in AA and pp collisions normalized on the number of binary nucleon- 
nucleon sub-collisions: 

dafjdpT , . 

TAA{b)a,^da;;/dpT ' ^ ^ 

where (7;^ = 42 mb is the inelastic non- diffract ive pp cross section at ^/s = 200 GeV. In the 
absence of medium-induced effects in the mid-rapidity region it should be Raa = 1 for high 
enough pt{^ 2 GeV/c). Such value of RAAsiml has been observed so for d+Au and peripheral 
Au-|-Au collisions, but not for central and semi-central Au-|-Au events, where Raa < 1 up to 
maximum measured transverse momenta Pt ~ 10 GeV/c. One can see from Figure 3 that our 
model calculations reproduce Pt- and centrality dependences of Raa [35] quite well. 

Another important tool to verify jet quenching is two-particle azimuthal correlation function 
C(A(/9) - the distribution over an azimuthal angle of high-p^ hadrons in the event with 2 
GeV/c < Pt < p%^^ relative to that for the hardest "trigger" particle with p!^'^ > 4 GeV/c. 
Figure 4 presents C(A(/9) in pp and in central Au+Au collisions (data from STAR [36]). Clear 
peaks in pp collisions at A(/9 = and A(/9 = tt indicate a typical dijct event topology. Note 
that almost the same pattern has been observed in d+Au and peripheral Au+Au collisions. 
However, for most central Au+Au collisions the peak near tt disappears. It can be interpreted 
as the observation of monojet events due to the "absorption" of one of the jets in a dense 
medium. Such event configuration corresponds to the situation when the dijct production 
vertex is close to the surface of the nuclear overlap region: then one partonic jet can escape the 
medium almost without re-interactions and then go to detectors, while second jet loses most of 
its initial energy due to a large number of rescatterings and therefore becomes unobscrvable [37]. 
Figure 4 demonstrates that measured suppression of azimuthal back-to-back correlations is well 
reproduced by our model (the same procedure of uncorrelated background subtraction as in [36] 
was applied). 

We leave beyond the scope of this paper the analysis of such important RHIC observables 
as the azimuthal anisotropy and particle ratios. Since these observables are very sensitive to 
the soft physics, in order to study them a more careful treatment of \ow-pT particle produc- 
tion than our simple approach is needed (the detailed description of space-time structure of 
freeze-out region, resonance decays, etc.) For example, our model can reproduce experimentally 
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measured [38, 39] py-dependence of the coefficient of azimuthal anisotropy V2 (the second har- 
monic of Fourier decomposition of particle azimuthal distribution) qualitatively, giving rapid 
hydrodynamical growth up to ~ 3 GeV/c with the subsequent saturation. However, the 
model calculations significantly underestimate the data at < 2 GeV/c. A solution of baryon- 
to-meson ratio "puzzle" [40, 41] is also beyond our consideration here. Further development of 
our model with special emphasis on the more detailed description of low-pr particle production 
is planed for the future. 

6 Conclusions 

The model of jet quenching in ultrarelativistic heavy ion collisions has been developed. It 
includes the generation of the hard parton production vertex according to the realistic nuclear 
geometry, rescattering-by-rescattering simulation of the parton path length in a dense matter, 
radiative and collisional energy loss per rescattering, final hadronization with the Lund string 
fragmentation model for hard partons and in-medium emitted gluons. The model is the fast 
Monte-Carlo tool implemented to modify a standard PYTHIA jet event. The model has been 
generalized to the case of the "full" heavy ion event (the superposition of soft, hydro- type state 
and hard multi-jets) using a simple and fast simulation procedure for soft particle production. 

The efficiency of the model is demonstrated basing on the numerical analysis of high-p^^ 
hadron production in Au+Au collisions at RHIC. The good fit of experimental data on r/- 
and Pt- spectra of hadrons for different event ccntralitics is achieved. The model is capable 
of reproducing main features of the jet quenching pattern at RHIC: the px dependence of the 
nuclear modification factor RaA) and the suppression of azimuthal back-to-back correlations. 
The further development of the model focusing on a more detailed description of low-pT particle 
production is planed for the future. 
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Figure 1 : The pseudorapidity distribution of charged hadrons in Au+ Au colhsions for three centrahty 
sets. The points are PHOBOS data [34], histograms are the model calculations. 
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Figure 2: The transverse momentum distribution of neutral pions in Au+Au collisions for three 
centrality sets. The points are PHENIX data [35], histograms are the model calculations. 
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Figure 3: The nuclear modification factor Raa (21) for neutral pions in Au+Au collisions for two 
centrality sets. The points are PHENIX data [35], histograms are the model calculations. 
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Figure 4: The azimuthal two-particle correlation function for pp and for central Au+Au collisions. 
The points are STAR data [35], dashed and solid histograms are the model calculations for pp and 
Au+Au events respectively. 
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